#### 0.准备好环境 ####
# Paper02-Sample_02_paper_成人非酒精性脂肪性肝病患者血清维生素 C 水平与肝纤维化的关系

# 引入包
library(haven)
library(plyr)
library(dplyr)
#加权分析用
library(survey)
#tableby
library(arsenal)
library(gtsummary)

# 4.3节用到包
library(tidyverse)
library(gtsummary)
library(tidyr) # drop_na 函数，用于快速去掉 NA

setwd("G:/BaiduNetdiskDownload/NHANES/")
#### 1. 定位数据模块和变量，获取源数据 ####
##### 1.1 DEMO-人口学数据提取 #####
demo.j <- read_xpt("2017-2018/Demographics/demo_j.xpt")
# colnames(demo.j)

### 2 提取研究所需要的变量
#  提取这些
# 年龄 RIDAGEYR
# 性别 RIAGENDR
# 种族 RIDRETH3
# 教育水平 DMDEDUC2
demo.data <- demo.j[,c('SEQN', 'RIDAGEYR', 'RIAGENDR', 'RIDRETH3', 'DMDEDUC2',"WTINT2YR", "SDMVPSU", "SDMVSTRA")]

##### 1.2 PAQ-运动数据提取 ##########
### 1.提取 Component文件
# https://wwwn.cdc.gov/Nchs/Nhanes/2017-2018/PAQ_J.htm
paq.j <- read_xpt("2017-2018/Questionnaire/paq_j.xpt")
### 2.提取研究所需要的变量
# 2023年2月22日14:54:40 这里不需要转换 直接拿取就可以了
# paq.data.file <- dplyr::bind_rows(list(paq.j))

paq.data <- paq.j[,c('SEQN', 'PAQ650','PAQ665')]


##### 1.3 BMI 数据提取 ##########
# https://wwwn.cdc.gov/Nchs/Nhanes/2017-2018/BMX_J.htm
### 1.提取 Component文件
bmx.j <- read_xpt("2017-2018/Examination/bmx_j.xpt") #注意是 Examination 的类别

### 2.提取研究所需要的变量 BMI-BMXBMI; 腰围-BMXWAIST
bmx.data <- bmx.j[,c('SEQN', 'BMXBMI')]

##### 1.4 SMQ-吸烟数据提取 ######
# https://wwwn.cdc.gov/Nchs/Nhanes/2017-2018/SMQ_J.htm
### 1.1 提取 Component文件
smq.j <- read_xpt("2017-2018/Questionnaire/smq_j.xpt")

### 2. 提取研究所需要的变量
# 是否吸烟至少100支-SMQ020; 现在是否吸烟-SMQ040; 多久前开始戒烟-SMQ050Q;
# 多久前开始戒烟的时间单位（天、周、月、年）-SMQ050U;
#smq.data.file <- dplyr::bind_rows(list(smq.g, smq.h))
# colnames(smq.j)
smq.data <- smq.j[,c('SEQN',"SMQ020", 'SMQ040', 'SMQ050Q')]

##### 1.5 糖尿病诊断  #####
# Participants with diabetes were identified as having any of the following: 
# (a) hemoglobin A1C concentration >= 6.5% ***or*** a fasting plasma glucose level >= 126 mg/dL [21]; 
# (b) for those who responded ‘yes’ to the question: ‘Doctor told you have diabetes?’ ***or*** ‘Taking insulin now?’.

ghb.j <- read_xpt("2017-2018/Laboratory/ghb_j.xpt")
ghb.data <- ghb.j[,c('SEQN','LBXGH')]

##### 1.5.2 空腹血糖（mg/dl）##### 
# https://wwwn.cdc.gov/Nchs/Nhanes/2017-2018/GLU_J.htm
glu.j <- read_xpt("2017-2018/Laboratory/glu_j.xpt")
# colnames(glu.j)
glu.data <- glu.j[,c('SEQN','LBXGLU')]

##### 1.5.3 是否有医生告知您患有糖尿病、现在是否使用胰岛素 ##### 
diq.j <- read_xpt("2017-2018/Questionnaire/diq_j.xpt")
# colnames(diq.j)
diq.data <- diq.j[,c('SEQN','DIQ010','DIQ050')]

#合并提取的四个维度的数据
diabetes <- plyr::join_all(list(ghb.data,glu.data,diq.data),by = "SEQN")
diabetes.a.index <- ifelse(diabetes$LBXGH>=6.5 | diabetes$LBXGLU >=126,1,NA)
diabetes.b.index <- ifelse(diabetes$DIQ010==1 | diabetes$DIQ050 ==1,1,NA)
# length(diabetes$LBXGH)
# length(diabetes$LBXGLU)
#这两个指标合并一起
diabetes.index <- ifelse(diabetes.a.index==1 | diabetes.b.index==1,1,0)
# 附件一起
diabetes$diabetes.index <- diabetes.index

# length(diabetes)

##### 1.6 丙氨酸转氨酶  ##########
# https://wwwn.cdc.gov/Nchs/Nhanes/2017-2018/BIOPRO_J.htm
alt.j <- read_xpt("2017-2018/Laboratory/biopro_j.xpt")
# colnames(alt.j)
alt.data <- alt.j[,c('SEQN',"LBXSATSI")]

##### 1.7 血清维生素C的含量 #####
# https://wwwn.cdc.gov/Nchs/Nhanes/2017-2018/VIC_J.htm
vic.j <- read_xpt("2017-2018/Laboratory/vic_j.xpt")
# colnames(vic.j)
vic.data <- vic.j[,c('SEQN',"LBXVIC","LBDVICSI")]

##### 1.8 高密度脂蛋白 #####
# https://wwwn.cdc.gov/Nchs/Nhanes/2017-2018/HDL_J.htm
hdl.j <- read_xpt("2017-2018/Laboratory/hdl_j.xpt")
# colnames(hdl.j)
hdl.data <- hdl.j[,c('SEQN',"LBDHDDSI","LBDHDD")]

##### 1.9 弹性成像检查  #####
# https://wwwn.cdc.gov/Nchs/Nhanes/2017-2018/LUX_J.htm
# 弹性成像检查状态-LUAXSTAT
# 可控衰减参数中位数-LUXCAPM
# 中位硬度-LUXSMED
lux.j <- read_xpt("2017-2018/Examination/lux_j.xpt")
# colnames(lux.j)
lux.data <- lux.j[,c('SEQN',"LUXCAPM","LUXSMED","LUAXSTAT")]

##### 1.10 饮酒数据提取 #####
# https://wwwn.cdc.gov/Nchs/Nhanes/2017-2018/DR1TOT_J.htm
dr1tot.j <- read_xpt('2017-2018/Dietary/dr1tot_j.xpt')
# colnames(dr1tot.j)
dr1tot.data <- dr1tot.j[,c('SEQN',"DR1TALCO")]

dr2tot.j <- read_xpt('2017-2018/Dietary/dr2tot_j.xpt')
dr2tot.j$DR2TALCO
# colnames(dr2tot.j)
dr2tot.data <- dr2tot.j[,c('SEQN',"DR2TALCO")]


alco.data <- merge(dr1tot.data, dr2tot.data)
# dim(alco.data)
# #View(alco.data)

##### 1.11 其他原因导致的肝脏疾病-viral hepatitis infection #####
# 乙型肝炎表面抗原
viral.hepbd.j <- read_xpt('2017-2018/Laboratory/hepbd_j.xpt')
viral.hepbd.data <- viral.hepbd.j[,c('SEQN',"LBDHBG")]

# 丙肝抗体\丙肝RNA
viral.hepc.j <- read_xpt('2017-2018/Laboratory/hepc_j.xpt')
viral.hepc.data <- viral.hepc.j[,c('SEQN',"LBXHCR", "LBDHCI")]

# 肝癌、自身免疫性肝炎
mcq.j <- read_xpt('2017-2018/Questionnaire/mcq_j.xpt')
liver.disease.data <- mcq.j[,c('SEQN',"MCQ510E", "MCQ230A", "MCQ230B", "MCQ230C")]


#### 2 提取分析相关变量（权重等，暂时为复现paper结果而提取） ####
##### 2.1 权重变量 ##### 
# 找到权重变量这块我是看不懂
# DEMO->MEC-Dietary, GLU(空腹血糖)
# 饮食中的权重：https://wwwn.cdc.gov/Nchs/Nhanes/2017-2018/DR1TOT_J.htm
# WTDRD1, WTDRD2

# GLU: https://wwwn.cdc.gov/Nchs/Nhanes/2017-2018/GLU_J.htm
# WTSAF2YR

# 权重的计算方式，需要按 Cycle 进行平均，这个仅有1个 Cycle，无需做处理
weight.data <- demo.j[,c('SEQN', 'WTMEC2YR')] 

##### 2.2 复杂抽样的其他变量-DEMO #####
survey.design.data <- demo.j[,c('SEQN', 'SDMVPSU', 'SDMVSTRA')]

#### 3.合并上述所有数据（把列拼接起来）-Output #### 
# 【4个# 用作第一层级】

##### 3.1 合并步骤1 & 2中提取的数据 #####
output <- plyr::join_all(list(demo.data, paq.data, bmx.data, smq.data,
                              ghb.data, glu.data, diq.data, diabetes,
                              alt.data, vic.data, hdl.data, lux.data, alco.data,
                              viral.hepbd.data, viral.hepc.data,liver.disease.data,
                              weight.data, survey.design.data), by='SEQN', type='full')

# #View(output) # 跑通到这一步，得到这篇文章所用到的全部原始数据
#### 4.根据文章的筛选策略进行筛选 ####
# 原始筛选策略：
#####  4.1缺少血清维生素 C 数据(n = 2514) #####
dim(output)

vc.exclude.index <- which(is.na(output$LBDVICSI))
length(vc.exclude.index)

#删除这些没有数据的 有- 是排除这些人没有- 是提取这些人
data.vc.exist <- output[-vc.exclude.index, ]
#dim(data.vc.exist) #6740

##### 4.2 排除 LUX 数据缺失以及不是 NAFLD 的人员 #####
# 弹性成像检查状态-LUAXSTAT
# 可控衰减参数中位数-LUXCAPM
# 中位硬度-LUXSMED
# table(output$LUAXSTAT)
# complete ineligible   not done    partial 
# 5494        258        156        493 
# 这块排除需要仔细阅读原文并且有点难啊这块什么意思
cap.exclude.index <- which(data.vc.exist$LUAXSTAT != 1 ## 弹性成像检查状态-LUAXSTAT
                           | is.na(data.vc.exist$LUXCAPM)  #可控衰减参数中位数-LUXCAPM
                           | is.na(data.vc.exist$LUXSMED)  # 中位硬度-LUXSMED
                           | data.vc.exist$LUXCAPM <248) # 可控衰减参数中位数-LUXCAPM
# length(cap.exclude.index) # 3945
data.vc.exist.cap.exist <- data.vc.exist[-cap.exclude.index,]
# dim(data.vc.exist.cap.exist) #[1] 2795
##### 4.3 排除其他原因导致的肝脏疾病 & 重要协变量的缺失 #####
# 乙型肝炎表面抗原-LBDHBG, # 丙肝抗体\丙肝RNA-LBDHCI, LBXHCR, # 肝癌MCQ230A-C、自身免疫性肝炎-MCQ510E
# 肝癌 MCQ230A-C
liver.cancer.index <- which(data.vc.exist.cap.exist$MCQ230A == 22|data.vc.exist.cap.exist$MCQ230B == 22|
                              data.vc.exist.cap.exist$MCQ230C == 22)

# 自身免疫性肝炎-Autoimmune hepatitis
autoimm.hepa.index <- which(data.vc.exist.cap.exist$MCQ510E == 5)

# 丙肝抗体\丙肝RNA-LBDHCI, LBXHCR
hepc.index <- which(data.vc.exist.cap.exist$LBDHCI == 1|data.vc.exist.cap.exist$LBXHCR == 1)

# 乙型肝炎表面抗原-LBDHBG
hepbd.index <- which(data.vc.exist.cap.exist$LBDHBG == 1)

# 过度饮酒
# 先计算2天的平均值，其中第二天为 NA 的，平均值也是 NA
total.alco <- apply(data.vc.exist.cap.exist[,c('DR1TALCO', 'DR2TALCO')], 1, mean)
data.vc.exist.cap.exist$total.alco <- total.alco
#View(data.vc.exist.cap.exist[,c('DR1TALCO', 'DR2TALCO', 'total.alco')])

# 把第二天为 NA 的值计算出来，用第一天的值作为平均值
day.2.na.index <- which(is.na(data.vc.exist.cap.exist$DR2TALCO))
# 看下替换前的结果
 #View(data.vc.exist.cap.exist[day.2.na.index,c('DR1TALCO', 'DR2TALCO', 'total.alco')])


data.vc.exist.cap.exist$total.alco[day.2.na.index] <- data.vc.exist.cap.exist$DR1TALCO[day.2.na.index]
# 看下替换后的结果
#View(data.vc.exist.cap.exist[day.2.na.index,c('DR1TALCO', 'DR2TALCO', 'total.alco')])

# 过度饮酒的定义：Excessive alcohol consumption 
excessive.alco.male <- ifelse(data.vc.exist.cap.exist$RIAGENDR == 1 & data.vc.exist.cap.exist$total.alco > 20, 1, 0)
excessive.alco.female <- ifelse(data.vc.exist.cap.exist$RIAGENDR == 2 & data.vc.exist.cap.exist$total.alco > 10, 1, 0)
table(excessive.alco.female)

excessive.alco.index <- which(excessive.alco.male == 1 | excessive.alco.female == 1, 1, 0)

# 汇总上述指标
other.cause.index <- c(autoimm.hepa.index, 
                       hepc.index, hepbd.index,
                       liver.cancer.index, 
                       excessive.alco.index)
# 最后结果的汇总
data.vc.exist.cap.exist.non.other.cause <- data.vc.exist.cap.exist[-other.cause.index,]
dim(data.vc.exist.cap.exist.non.other.cause) #2374  

# covariate 上的缺失数据 也就是协变量缺失的人
paper.data <- subset.data.frame(data.vc.exist.cap.exist.non.other.cause,
                                (!is.na(LBDVICSI)) & # 维C，0，之前排除过缺失
                                  (!is.na(LBXSATSI)) & # 丙氨酸转氨酶, 28
                                  (!is.na(LBDHDDSI)) & # 高密度脂蛋白, 15
                                  (!is.na(LUXCAPM)) & # 可控衰减参数中位数-LUXCAPM，0，之前排除过
                                  (!is.na(LBXGH)) & # 糖化血红蛋白A1C(%), 2
                                  (!is.na(total.alco)) & # 24小时回忆的饮酒量
                                  (!is.na(RIDAGEYR)) & # 年龄
                                  (!is.na(RIAGENDR)) & # 性别
                                  (!is.na(RIDRETH3)) & #种族
                                  (!is.na(DMDEDUC2))  # 教育程度
)


dim(paper.data) # [1] 1926   36
# 直接使用现成数据 过滤太麻烦了
#write.csv(paper.data,'2.9.csv')
#A median liver stiffness of 8.2 kPa was used to define cases of significant fibrosis (F2) in this study 
# 用于区分定义严重纤维化 分值
paper.data$fibrosis_group = ifelse(paper.data$LUXSMED > 8.2 , "YES" , "No")
dim(paper.data[which(paper.data$fibrosis_group=='Yes'),])
dim(paper.data[which(paper.data$fibrosis_group!='Yes'),])


paper.data$Sex <- ifelse(paper.data$RIAGENDR==1,'male','female')
paper.data$Raceethnicity <- recode_factor(paper.data$RIDRETH3,
                                                                `1` = 'Mexican American',
                                                                `2` = 'Other Hispanic',
                                                                `3` = 'Non-Hispanic White',
                                                               `4` = 'Non-Hispanic Black',
                                                                `6` = 'Other/multiracial',
                                                              `7` = 'Other/multiracial')

paper.data$Educationlevel <- ifelse(paper.data$DMDEDUC2<=3,'<=High school','> High school')

paper.data$BMIgroup <- ifelse(paper.data$BMXBMI < 25, '<25',
                               ifelse(paper.data$BMXBMI >=25 & paper.data$BMXBMI < 30, '25-30',
                                      ifelse(paper.data$BMXBMI >=30, '>=30',
                                             'NA')))

paper.data$rpa <- ifelse(paper.data$PAQ650 ==1|paper.data$PAQ665==1, 'YES','No')

#吸烟原文判断 如果
# Participants were asked whether they had ever smoked 100 cigarettes in their lifetime 
# and whether they smoked currently to identify current
# and former smokers. Participants were defined as former
# smokers if they did not smoke currently but had ever
# smoked  100 cigarettes in the past
# 受访者问到一生中是不是吸过100 只烟 以及当前是不是吸烟来判断现在吸烟和以前吸烟 下面是解释
# 被定义为曾经吸烟者的情况是现在不吸烟但是在过去曾经吸过100只烟

# 一点也不吸烟 没说怎么判断 用我的方法得出的数量与论文不一致那么就是去除下面两种情况的都是从不吸烟
dim(paper.data[which(paper.data$SMQ040==3),])
# 当前吸烟
dim(paper.data[which(paper.data$SMQ040==1|paper.data$SMQ040==2),])
#曾经吸烟
dim(paper.data[which(paper.data$SMQ040==3 ),])


paper.data$smoke_group = ifelse((paper.data$SMQ040==1|paper.data$SMQ040==2) , "Current",
                                ifelse((paper.data$SMQ040==3) , "Former" ,'Never'))
paper.data$smoke_group[which(is.na(paper.data$smoke_group))] <-'Never'

dim(paper.data[which(paper.data$smoke_group=='Never' ),])


paper.data$diabetesindex <-ifelse(paper.data$diabetes.index==1,'YES')
paper.data$diabetesindex[which(is.na(paper.data$diabetes.index))] <-'No'


paper.data$fibrosis_groupint <- ifelse(paper.data$fibrosis_group=='YES',1,0)

paper.data$LBXVIC.quantile.var <- cut(paper.data$LBXVIC,
                                      breaks = quantile(paper.data$LBXVIC),
                                      labels = c('Q1', 'Q2', 'Q3', 'Q4'))



groupdata <- paper.data %>% select(RIDAGEYR,Sex,Raceethnicity,Educationlevel,BMIgroup,rpa,smoke_group,diabetesindex,LUXSMED,
                                   LBXVIC,LBXSATSI,LBDHDD,LUXCAPM,diabetes.index,
                                   fibrosis_group,SDMVPSU,SDMVSTRA,WTINT2YR,fibrosis_groupint,LBXVIC.quantile.var)
#write.csv(paper.data,'2.9csv.csv')
#上面不需要 与下面的复杂抽样没什么区别 可以直接放入原始数据

# 加权线性回归方程测试P值

NHANES_design <- svydesign(
  data = paper.data,
  ids = ~SDMVPSU,
  strata = ~SDMVSTRA,
  nest = TRUE,
  weights = ~WTINT2YR,survey.lonely.psu = "adjust"
)


tbl_svysummary(NHANES_design,
              missing = 'no', 
              # 通过属性分组
              by = fibrosis_group,
              # 分区yes 需要YES 和NO 如果采用Yes 和No 不知道是不是有变量冲突导致分类不显示 只显示全部的
              # 进行设置显示维度
              include = c(RIDAGEYR,Sex,Raceethnicity,Educationlevel,BMIgroup,rpa,smoke_group,diabetesindex,
                          LBXVIC,LBXSATSI,LBDHDD,LUXCAPM),
              # 批量给变量起别名
              label = list(RIDAGEYR ~ "Age (years)",
                           LUXCAPM ~ "Median CAP (dB/m)",
                           Sex ~ "Gender, n (%)a",
                           rpa ~ "Recreational physical activity, n (%)a",
                           Raceethnicity ~ "Race/ethnicity, n (%)a",
                           Educationlevel ~ "Education level, n (%)a",
                           BMIgroup ~ "BMI group, n (%)a",
                           smoke_group ~ "Smoking status, n (%)a",
                           LBXVIC ~ "Serum vitamin C (mg/dL)",
                           diabetesindex ~ "Diabetes, n (%)a",
                           LBXSATSI ~ "ALT (IU/L)",
                           LBDHDD ~ "HDL-cholesterol (mg/dL)"),
              statistic = list(
                # 分别对应数值型和分类变量 后是需要显示表达式
                all_continuous() ~ "{mean}±{sd}",
                # 样本是加权倒置n也是加权的 下面参数是显示
                all_categorical() ~ "{n_unweighted} ({p}%)"
              ),
              digits = list(all_continuous() ~ 2
                            # ,all_categorical() ~ 3
                            )
              )%>%
   add_n(
     statistic = "{N_nonmiss_unweighted}",
     col_label = "**N**"
   
   ) %>% # 添加非NA观测值个数
  # add_p() %>%
  # 这里p值不能用p 根据原文说的是加权得来的 需要用 tbl_svysummay进行 这里除了p值 其他是正确的
  add_p() %>% # 添加P值
  add_overall() %>%
  #modify_header(all_stat_cols() ~ "**level**, N = {n_unweighted} ({stype_percent(p)}%)" ) %>%
  # 由原来的加权显示改成非加权显示原始值 
  modify_header(stat_1 ~ "**No** (n={n_unweighted})") %>%
  modify_header(stat_2 ~ "**Yes** (n={n_unweighted})") %>%
  modify_spanning_header(c("stat_1", "stat_2") ~ "**Significant fibrosis**")%>%as_flex_table()%>%flextable::save_as_html(path = '2.9tab1.html')
  # 编辑表头测试
  #modify_spanning_header(stat_1 ~ NA, update = all_stat_cols() ~ "**Significant fibrosis**") %>%
  #modify_footnote(update = all_stat_cols() ~ "测试编辑表格下面值")
  #%>%
  #bold_labels()  #label 粗体


  
# show_header_names(svytable1)
# 
# table(groupdata$fibrosis_group)
# # 
#  
# svytable1 
  
#成人非酒精性脂肪性肝病患者血清维生素 C 水平与肝纤维化的关系 y变量是肝纤维化 x 变量是血清维生素C 水平

# tbl_summary(groupdata,
#             # 通过属性分组
#             by = fibrosis_group,
#             # 分区yes 需要YES 和NO 如果采用Yes 和No 不知道是不是有变量冲突导致分类不显示 只显示全部的
# 
#             # 进行设置显示维度
#             include = c(RIDAGEYR,Sex,Raceethnicity,Educationlevel,BMIgroup,rpa,smoke_group,diabetesindex,
#                                                       LBXVIC,LBXSATSI,LBDHDD,LUXCAPM),
#             # 批量给变量起别名
#             label = list(RIDAGEYR ~ "Age (years)",
#                          LUXCAPM ~ "Median CAP (dB/m)",
#                          Sex ~ "Gender, n (%)a",
#                          rpa ~ "Recreational physical activity, n (%)a",
#                          Raceethnicity ~ "Race/ethnicity, n (%)a",
#                          Educationlevel ~ "Education level, n (%)a",
#                          BMIgroup ~ "BMI group, n (%)a",
#                          smoke_group ~ "Smoking status, n (%)a",
#                          LBXVIC ~ "Serum vitamin C (mg/dL)",
#                          diabetesindex ~ "Diabetes, n (%)a",
#                          LBXSATSI ~ "ALT (IU/L)",
#                          LBDHDD ~ "HDL-cholesterol (mg/dL)"),
# 
#             statistic = list(
#               # 分别对应数值型和分类变量 后是需要显示表达式
#               all_continuous() ~ "{mean}±{sd}",
#               all_categorical() ~ "{n} ({p}%)"
#             ),
#             digits = all_continuous() ~ 3
#             )%>%
#   add_n() %>% # 添加非NA观测值个数
#   # 这里p值不能用p 根据原文说的是加权得来的 需要用 tbl_svysummay进行 这里除了p值 其他是正确的
#   add_p(
#     # https://www.danieldsjoberg.com/gtsummary/reference/tests.html
#     list(all_continuous() ~ "t.test", all_categorical() ~ "chisq.test")
#     ) %>% # 添加P值
#   add_overall() %>%
#   modify_spanning_header(c("stat_1", "stat_2") ~ "**Significant fibrosis**") %>%
#   #modify_header(label = "**Variable**") %>% # 标签列header
#   bold_labels()  #label 粗体
# 现在p值有问题 我应该是算原始值 而不是转换后,这只是我猜想



#### table2 进行测试####
# 线性回归模型进行测试
# Serum vitamin C
# NHANES_designfg1 <-subset(NHANES_design,fibrosis_groupint==1)
# model1
glmsvc <- svyglm (fibrosis_groupint ~ LBXVIC+LBXVIC.quantile.var ,design = NHANES_design,family=quasibinomial )
model1 <-tbl_regression(glmsvc,exponentiate = T,label = list(LBXVIC ~ 'Serum vitamin C (mg/dL)',LBXVIC.quantile.var ~ 'Serum vitamin C (mg/dL, quartile)') )%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")
# model2

#  模型2 年龄性别和种族
# RIDAGEYR +  Sex + Raceethnicity

# 多元由原来的分类变量 试试连续型变量
glmsvcmodel2 <- svyglm (fibrosis_groupint ~ LBXVIC+LBXVIC.quantile.var+RIDAGEYR +  Sex + Raceethnicity ,design =NHANES_design,family=quasibinomial )

model2<-tbl_regression(glmsvcmodel2,include = c(LBXVIC,LBXVIC.quantile.var),label = list(LBXVIC ~ 'Serum vitamin C (mg/dL)',LBXVIC.quantile.var ~ 'Serum vitamin C (mg/dL, quartile)') ,exponentiate = T)%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")


# model3 
#模型3  年龄、性别、种族/族裔、BMI、教育水平、吸烟状况、休闲体育活动、高密度脂蛋白胆固醇和糖尿病
# RIDAGEYR +  Sex + Raceethnicity + BMIgroup + Educationlevel+ smoke_group+ rpa+  LBDHDD+  diabetes.index

glmsvcmodel3 <- svyglm(fibrosis_groupint ~ LBXVIC+LBXVIC.quantile.var+RIDAGEYR +
                          Sex + Raceethnicity + BMIgroup + Educationlevel +
                          smoke_group+ rpa+  LBDHDD + diabetes.index ,design = NHANES_design,family=quasibinomial)

model3<-tbl_regression(glmsvcmodel3,include = c(LBXVIC,LBXVIC.quantile.var),label = list(LBXVIC ~ 'Serum vitamin C (mg/dL)',LBXVIC.quantile.var ~ 'Serum vitamin C (mg/dL, quartile)') ,exponentiate = T)%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")


alllei2<-tbl_merge(tbls = list(model1,model2,model3),tab_spanner = c('Model 1OR (95% CI)','Model 2 OR (95% CI)','Model 3OR (95% CI)'))
alllei2%>%as_flex_table()%>%flextable::save_as_html(path = "2.9table2.html")
#### table3  进行测试####
#model1 
# model2
# RIDAGEYR +  Sex + Raceethnicity
# model3 
#模型3  年龄、性别、种族/族裔、BMI、教育水平、吸烟状况、休闲体育活动、高密度脂蛋白胆固醇和糖尿病
# RIDAGEYR +  Sex + Raceethnicity + BMIgroup + Educationlevel+ smoke_group+ rpa+  LBDHDD+  diabetes.index
#model1  这里面有个问题 如果你设置的截取数据在你分析协变量的维度里面是不能是已经拆分出固定的值 就比如数据里面性别选择了男 那么协变量不能有性别维度进行分析会报错

glmsvcmodel3male <- svyglm(fibrosis_groupint ~ LBXVIC
                           ,design = subset(NHANES_design,Sex =='male'))%>% 
  tbl_regression(exponentiate = TRUE, include = c(LBXVIC),label = list(LBXVIC ~ 'male'))%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")


glmsvcmodelfemale <- svyglm(fibrosis_groupint ~ LBXVIC,
                            design = subset(NHANES_design,Sex =='female')) %>% 
  tbl_regression(exponentiate = TRUE, include = c(LBXVIC),
                 label = list(LBXVIC ~ 'female'))%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")
# Overweight

glmsvcmodel3obese <- svyglm(fibrosis_groupint ~ LBXVIC
                            ,design = subset(NHANES_design,BMIgroup =='>=30'))%>% 
  tbl_regression(exponentiate = TRUE, include = c(LBXVIC),label = list(LBXVIC ~ 'obese'))%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")

glmsvcmodel3underweight <- svyglm(fibrosis_groupint ~ LBXVIC
                                  ,design = subset(NHANES_design,BMIgroup =='<25'))%>% 
  tbl_regression(exponentiate = TRUE, include = c(LBXVIC),label = list(LBXVIC ~ 'Under/normal weight'))%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")

glmsvcmodel3Overweight <- svyglm(fibrosis_groupint ~ LBXVIC
                                 ,design = subset(NHANES_design,BMIgroup =='25-30'))%>% 
  tbl_regression(exponentiate = TRUE, include = c(LBXVIC),label = list(LBXVIC ~ 'overweight'))%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")

model1<- tbl_stack(tbls = list(glmsvcmodel3male,glmsvcmodelfemale,glmsvcmodel3underweight,glmsvcmodel3Overweight,glmsvcmodel3obese),group_header = c('Stratified by gender','Stratified by gender','Stratified by BMI','Stratified by BMI','Stratified by BMI'))

#model2 
# RIDAGEYR +  Sex + Raceethnicity

glmsvcmodel3male <- svyglm(fibrosis_groupint ~ LBXVIC+RIDAGEYR + Raceethnicity
                           ,design = subset(NHANES_design,Sex =='male'))%>% 
  tbl_regression(exponentiate = TRUE, include = c(LBXVIC),label = list(LBXVIC ~ 'male'))%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")


glmsvcmodelfemale <- svyglm(fibrosis_groupint ~ LBXVIC+RIDAGEYR  + Raceethnicity,
                            design = subset(NHANES_design,Sex =='female')) %>% 
  tbl_regression(exponentiate = TRUE, include = c(LBXVIC),
                 label = list(LBXVIC ~ 'female'))%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")
# Overweight

glmsvcmodel3obese <- svyglm(fibrosis_groupint ~ LBXVIC+RIDAGEYR +  Sex + Raceethnicity
                            ,design = subset(NHANES_design,BMIgroup =='>=30'))%>% 
  tbl_regression(exponentiate = TRUE, include = c(LBXVIC),label = list(LBXVIC ~ 'obese'))%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")

glmsvcmodel3underweight <- svyglm(fibrosis_groupint ~ LBXVIC+RIDAGEYR +  Sex + Raceethnicity
                                  ,design = subset(NHANES_design,BMIgroup =='<25'))%>% 
  tbl_regression(exponentiate = TRUE, include = c(LBXVIC),label = list(LBXVIC ~ 'Under/normal weight'))%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")

glmsvcmodel3Overweight <- svyglm(fibrosis_groupint ~ LBXVIC+RIDAGEYR +  Sex + Raceethnicity
                                 ,design = subset(NHANES_design,BMIgroup =='25-30'))%>% 
  tbl_regression(exponentiate = TRUE, include = c(LBXVIC),label = list(LBXVIC ~ 'overweight'))%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")



model2<- tbl_stack(tbls = list(glmsvcmodel3male,glmsvcmodelfemale,glmsvcmodel3underweight,glmsvcmodel3Overweight,glmsvcmodel3obese),group_header = c('Stratified by gender','Stratified by gender','Stratified by BMI','Stratified by BMI','Stratified by BMI'))


# model3 
# RIDAGEYR +  Sex + Raceethnicity + BMIgroup + Educationlevel+ smoke_group+ rpa+  LBDHDD+  diabetes.index
glmsvcmodel3male <- svyglm(fibrosis_groupint ~ LBXVIC+RIDAGEYR  + Raceethnicity + BMIgroup + Educationlevel+ smoke_group+ rpa+  LBDHDD+  diabetes.index
                           ,design = subset(NHANES_design,Sex =='male'))%>% 
  tbl_regression(exponentiate = TRUE, include = c(LBXVIC),label = list(LBXVIC ~ 'male'))%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")


glmsvcmodelfemale <- svyglm(fibrosis_groupint ~ LBXVIC+RIDAGEYR  + Raceethnicity + BMIgroup + Educationlevel+ smoke_group+ rpa+  LBDHDD+  diabetes.index,
                            design = subset(NHANES_design,Sex =='female')) %>% 
  tbl_regression(exponentiate = TRUE, include = c(LBXVIC),
                 label = list(LBXVIC ~ 'female'))%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")
# Overweight

glmsvcmodel3obese <- svyglm(fibrosis_groupint ~ LBXVIC+RIDAGEYR +  Sex + Raceethnicity  + Educationlevel+ smoke_group+ rpa+  LBDHDD+  diabetes.index
                            ,design = subset(NHANES_design,BMIgroup =='>=30'))%>% 
  tbl_regression(exponentiate = TRUE, include = c(LBXVIC),label = list(LBXVIC ~ 'obese'))%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")

glmsvcmodel3underweight <- svyglm(fibrosis_groupint ~ LBXVIC+RIDAGEYR +  Sex + Raceethnicity  + Educationlevel+ smoke_group+ rpa+  LBDHDD+  diabetes.index
                                  ,design = subset(NHANES_design,BMIgroup =='<25'))%>% 
  tbl_regression(exponentiate = TRUE, include = c(LBXVIC),label = list(LBXVIC ~ 'Under/normal weight'))%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")

glmsvcmodel3Overweight <- svyglm(fibrosis_groupint ~ LBXVIC+RIDAGEYR +  Sex + Raceethnicity  + Educationlevel+ smoke_group+ rpa+  LBDHDD+  diabetes.index
                                 ,design = subset(NHANES_design,BMIgroup =='25-30'))%>% 
  tbl_regression(exponentiate = TRUE, include = c(LBXVIC),label = list(LBXVIC ~ 'overweight'))%>%modify_column_hide(p.value)%>%modify_header(estimate ~ 'exp(Beta) (95% CI)')%>%modify_column_merge(pattern = "{estimate} ({ci})")


model3<- tbl_stack(tbls = list(glmsvcmodel3male,glmsvcmodelfemale,glmsvcmodel3underweight,glmsvcmodel3Overweight,glmsvcmodel3obese),group_header = c('Stratified by gender','Stratified by gender','Stratified by BMI','Stratified by BMI','Stratified by BMI'))

alllei3<-tbl_merge(tbls = list(model1,model2,model3),tab_spanner = c('Model 1OR (95% CI)','Model 2 OR (95% CI)','Model 3OR (95% CI)'))
alllei3%>%as_flex_table()%>%flextable::save_as_html(path = "2.9table3.html")


















